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DIGITAL  IMAGE  PROCESSING  BY  USE  OF 
LOCAL  STATISTICS 


INTRODUCTION 

Image  processing  on  digital  computers  has  been  gaining  in  acceptance  in  recent  years 
[1-3] . Early  techniques  in  image  processing  concentrated  mostly  on  procedures  that  were 
carried  out  computationally  in  the  frequency  domain  (Fourier  or  Walsh),  which  was  a 
natural  extension  of  one-dimensional  linear  signal  processing  theory.  In  due  course  it  be- 
came well  known  that  computing  a two-dimensional  transform  for  a large  data  array  is  a 
very  time-consuming  activity  even  with  fast  transform  techniques  on  large  computers. 

Hence,  implementation  of  frequency-domain  procedures  for  real-time  processing  of  images 
appears  less  than  promising.  More  recent  works  based  on  an  application  of  the  Kalman  filter- 
ing algorithm  [4]  or  Bayesian  estimation  extended  to  two-dimensional  arrays  led  to  the  con- 
cept of  a recursive  ..itering  algorithm  [5,6] . The  power  of  recursive  algorithms  for  real-time, 
one-dimensional  signal  processing  is  well  established.  However,  when  applied  to  a two- 
dimensional  array,  the  algorithm  operates  in  the  spatial  domain  in  which  pixels  have  to  be 
processed  sequentially.  As  a consequence,  the  procedure  is  no  longer  computationally  ef- 
ficient and  loses  its  attraction  for  real-time  processing. 

Algorithms  developed  in  this  report  share  a particular  characteristic  in  that  each  pixel 
can  be  processed  separately  without  waiting  for  its  neighboring  pixels  to  be  processed.  This 
characteristic  permits  a direct  implementation  of  these  edgorithms  for  real-time  image  proc- 
essing. Applying  local  statistics  to  image  processing  is  not  a new  idea.  Ketcham,  et  al.  [7] 
used  the  entire  local  histogram  for  real-time  image  enhancement,  and  Wallis  [8]  applied 
local  mean  and  variance  to  filter  out  scan-line  noise,  with  striking  results.  This  paper  extends 
this  family  of  algorithms  to  contrast  enhancement,  and  noise  filtering.  Both  additive  white 
noise  and  multiplicative  white  noise  cases  are  considered.  Most  additive  noise  filtering 
approaches  utilize  fast  Fourier  transform,  convolution  or  recursive  algorithms.  In  the  trans- 
form and  convolution  methods,  the  autocorrelation  between  pixels  is  either  assumed  or 
approximated,  and  in  the  recursive  algorithm,  a linear  causal  or  semicausal  autoregressive 
image  model  is  generally  assumed.  The  techniques  based  on  the  use  of  local  mean  and  vari- 
ance described  in  this  report  deviate  from  these  approaches.  The  basic  assumption  is  that 
the  sample  mean  and  variance  of  a pixel  is  equal  to  the  local  mean  and  variance  of  all  pixels 
within  a fixed  range  surrounding  it.  The  validity  of  this  assumption  is  debatable,  but  so  are 
most  other  statistical  image  representations  encountered  in  the  current  practice.  In  the  ad- 
ditive noise-filtering  case,  the  a priori  mean  (variance)  of  the  estimated  image  is  calculated 
as  the  difference  between  the  mean  (variance)  of  the  noise-corrupted  image  and  the  mean 
(variance)  of  the  noise  by  itself.  This  technique  is  extended  to  include  multiplicative  noise 
filtering  and  also  the  case  involving  both  multiplicative  and  additive  noise.  Although  this 
simple  approach  may  not  have  the  mathematical  elegance  and  sophistication  of  a few  other 
techniques,  our  experimental  results  and  those  reported  by  Wallis  [8]  indicate  that  the  local 
mean  and  variance  technique  is  a very  effective  tool  in  both  contrast-stretching  and  noise- 
filtering applications. 
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Let  Xi  j be  the  brightness  of  the  pixel  (i,  j)  in  a two-dimensional  N X N image.  The  local 
mean  and  variance  are  calculated  over  a {2n  + l)(2m  + 1)  window.  The  local  mean  is  defined 
as 


1 


n + i 


m+j 


(2n  + l)(2m  + 1) 


E E 


k^n-V.  « = 


m -j 


(1) 


Similarly,  the  local  variance  is 


m +j 
v-i 


(2n  + l)(2m  + l) 


n +1 

E ^ 

!=n-i  V.  = m -j 


m^jy 


(2) 


In  this  report,  the  next  section  is  devoted  to  spatial  contrast  enhancement  for  which 
only  the  local  mean  is  required.  The  third  section  deals  with  additive  noise  filtering;  the 
fourth  section  treats  images  corrupted  by  multiplicative  noise;  and  the  following  section 
extends  the  results  of  the  third  emd  fourth  sections  to  images  corrupted  by  both  additive 
and  multiplicative  noise.  In  the  fifth  section,  a simplified  algorithm  is  discussed  to  approxi- 
mate the  local  mean  and  variance  and  future  research  using  local  statistics  is  briefly  men- 
tioned. 


SPATIAL  CONTRAST  ENHANCEMENT 

Grey -level  rescaling,  high-pass  filtering,  and  histogram  redistribution  [9]  are  the  most 
popular  techniques  in  image  contrast  enhancement.  Wallis  [8]  suggested  an  algorithm  based 
on  local  mean  and  variance  in  which  each  pixel  is  required  to  have  a “desirable”  local  mean 
and  a “desirable”  local  variance  , such  that 


where  in  Eq.  (3),  m^  j and  u,  y are  local  mean  and  variance  as  given  by  Eqs.  (1)  and  (2).  It  is 
easy  to  verify  that  the  x'j  j h^  a mean  and  variance  u^,  if  we  consider  m,  y and  u,  y as 
the  true  mean  and  variance  of  x,  y.  The  main  drawback  of  this  technique  is  that  it  tends  to 
enhance  subtle  details  at  the  expense  of  the  principal  features,  which  are  lost  in  the  process. 
Figure  1 shows  the  original  image  and  the  image  processed  by  Eq.  (3).  The  river  in  the 
original  image  and  other  large  objects  are  difficult  to  recognize  in  the  processed  image.  To 
compensate  for  this,  we  designed  an  algorithm  such  that  a pixel  x,  y will  maintain  its  local 
mean,  and  yet  permits  its  variance  to  be  modified  by  a factor  of  its  originetl  loceil  variance. 
The  modified  algorithm  is 


^ij  " ^i,J  * (4) 

where  k,  the  gain,  is  the  ratio  of  new  local  standard  deviation  to  the  original  standard  devia- 
tion. This  scheme  heis  a distinct  computational  adveintage  over  the  scheme  rooted  in  Eq.  (3). 
The  computation  of  local  variance  u,  y is  not  required  and  only  m,  y need  be  computed.  If 
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fe  > 1,  the  image  will  be  sharpened  as  if  acted  upon  by  a high-pass  filter.  If  0 < fe  < 1,  the 
image  will  be  smoothed  as  if  passed  through  a low-pass  filter.  In  the  extreme  case,  k = 0 and 
X-  j is  equal  to  its  local  mean  y. 

This  algorithm  can  be  easily  extended  to  enhance  images  with  an  ill-distributed  histo- 
gram. Let^(x)  be  the  grey-level  rescaling  function  [9] ; then  Eq.  (4)  is  modified  to 

x'l  j = gim/  j)  + k(Xjj  -nijj).  (5) 

In  the  case  of  a linear  grey -level  stretch,  g(x ) is  a linear  relation.  Several  images  are  processed 
using  Eq.  (5)  and  are  shown  in  the  remainder  of  Fig.  1.  The  linear  function  g(x},  used  in  this 
picture,  yields  an  effective  contrast  stretch  in  both  the  highlights  and  the  dark  areas  of  the 
image.  The  window  size  of  3 X 3 or  5 X 5 is  large  enough  to  carry  out  contrast  enhance- 
ment. For  noise  filtering  (to  be  discussed  in  later  sections),  however,  a 7 X 7 or  larger 
dimensional  window  is  more  desirable.  All  images  of  Fig.  1 are  processed  with  a 5 X 5 
window.  Figure  Ic  shows  that  for  fe  = 2,  both  the  gross  features  and  subtle  details  are 
enhanced  in  the  same  proportion.  The  case  fe  = 1 has  no  effect  on  the  image.  For  k = 0.5, 
the  image  is  blurred  as  if  passed  through  a low-pass  filter.  For  fe  = 0,  the  image  is  averaged 
over  the  5X5  window. 


ADDITIVE  NOISE  FILTERING 

In  this  section,  attention  is  focused  on  filtering  of  images  degraded  by  white  additive 
noise.  Most  current  approaches  to  this  problem  employ  frequency-domain  techniques, 
direct  inversion,  or  recursive  Kalman  filtering.  In  this  report,  we  derive  a very  simple 
algorithm  based  on  the  use  of  local  mean  and  variance  of  an  image.  Let  z,  y be  the  degraded 
pixel  X,- j:  then 


where  tu,  y is  a white  random  sequence,  with 

£[w,j]  = 0 

and 

where  6,  is  the  Kronecker  delta  function  and  E is  the  expectation  operator.  (The  distribu- 
tion of  lu,  y is  irrelevant  in  this  derivation.)  The  problem  is  to  estimate  x,-  y,  given  z,-  y and 
the  noise  statistics. 

From  Eq.  (6),  we  have 


and 


(7) 

(8) 
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In  most  filtering  techniques,  the  a priori  mean  and  variance  of  Xj  j are  derived  from  an 
assumed  autocorrelation  model  or,  recursively,  from  an  autoregressive  Markov  sequence.  In 
either  method,  it  is  an  approximation.  Assume  thatTc,-  y and  Q,  y are  the  a priori  mean  and 
variance  of  x,  y,  which  in  turn  are  approximated  by  the  local  mean  and  variance  of  Eqs.  (7) 
and  (8).  Under  this  assumption,  it  is  very  easy  to  obtain  a filtering  algorithm  either  by 
minimizing  the  mean-square  error  or  by  weighted  least-square  estimation  [10] . Both 
methods  will  produce  the  same  algorithm.  The  estimated  x,-  y,  x,  is  computed  by 


where  the  gain  is  given  by 


= 


QiJ 


+ ai 


(10) 


Equation  (9)  can  be  written  as 

Xj  j - (l-ki  j)Xj  j + kj  j{Zj  j).  (11) 

Since  Q,  y and  Oj  are  both  positive,  fe,-  y will  lie  between  0 and  1.  A simple  intuitive  interpre- 
tation is  that,  for  a low  signal-to-noise  ratio  region,  Q,-  y is  small  compared  with  o^,  fe,  y ^ 0, 
and  the  estimated  Xj  y is'x;  . Conversely,  for  ahi^  si^al-to-noise  ratio  region,  Q,-  y is  much 
larger  than  Oj , fe,-  y ^ 1 , and  x,-  y *=  z,  y,  the  corrupted  pixel.  The  use  of  different  window 
sizes  will  greatly  affect  the  quality  o^  processed  images.  If  the  window  is  too  small,  the 
noise-filtering  algorithm  is  not  effective.  If  the  window  is  too  large,  subtle  details  of  the 
image  will  be  lost  in  the  filtering  process.  Our  experiments  indicate  that  a 7 X 7 window  is 
a fairly  good  choice.  All  images  presented  in  this  and  later  sections  are  processed  by  the 
7X7  window. 

Figure  2 shows  the  original  image,  the  image  contaminated  with  additive  noise  and  the 
estimated  image  produced  by  the  local  mean  and  variance  algorithm.  Clearly,  in  a smooth 
area,  the  pixel  is  averaged  over  the  window  and  in  a high-contrast  or  edge  area,  the  noise- 
corrupted  pixels  are  weighted  higher  than  their  local  mean  value.  Figures  2d,  2e  and  2f  are 
the  plots  of  intensity  along  a specific  scan  line  for  the  original,  the  noise-corrupted,  and  the 
processed  images,  respectively.  This  algorithm  works  equally  well  for  an  image  corrupted  by 
a Gaussian  noise.  Results  for  the  latter  case  eure  shown  in  Figs.  3a  and  3b. 


MULTIPLICATIVE  NOISE  FILTERING 

Images  containing  multiplicative  noise  have  the  characteristic  that  the  brighter  the 
area  the  noisier  it  is.  Mathematically,  the  degraded  pixel  can  be  represented  by 

^iJ  = • (12) 
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(a)  Original 


(b)  Original  plus  additive.  Uniform 
noise  (-30,  30) 


Ic)  Additive  noise  removed  with  (7  « 7)  mesh 
O^.  « 300 


<d)  Original  intensity  profile  along 
a scan  line 


Figure  2 
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Figure  3 


where 

= “/J 

EUUjj  -w'fc.t)]  = 

Nahi  and  Naraghi  [11]  treat  this  problem  via  the  Kalman-Bucy  approach,  which  necessitates 
solving  the  nonlinear  estimation  problem  by  numerical  integration.  In  this  report  the  non- 
linearity in  Eq.  (12)  is  treated  differently.  An  optimal  linear  approximation  of  Eq.  (12)  is 
used  to  produce  a filtering  algorithm  similar  to  that  for  the  additive  noise  case.  Experi- 
mental results  show  that  the  derived  algorithm  is  a very  promising  one. 

Let 


z'lj  = AXj  j + BUi  j + C,  (13) 

where  A,  B,  and  C are  not  random  variables.  They  are  to  be  chosen  to  minimize  the  mean- 
square  error  between  Zj  j and  z-  j and  also  to  make  z'j  j an  unbia.sed  estimate  of  Zj  j.  For 
z\  j to  be  the  unbiased  estimate  of  Zj  j,  we  must  have 

or 

C = j - AYj  j - BUi  j.  (14) 

Substituting  Eq.  (14)  into  (13)  and  forming  the  mean-square  error,  we  arrive  at  the  per- 
formance index  to  be  minimized. 
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J £[A(Xj- y *i,y^  ^ ~ ~ ' 

Upon  carrying  out  the  necessary  mathematical  procedures,  we  obtain  the  relation 

It  is  not  surjirising  to  find  that  Eq.  (15)  actually  is  the  first-order  Taylor  series  expansion  of 
2,  j about 

The  a priori  mean  and  variance  of  a:,  j are  computed  from  Eq.  (12)  and  are  given  by 


Var(2,.  ,)  + zfj 


in  which  Var(2,  y)  is  the  variance  ofzi  j.  The  quantities T,  j and  Var(2,  y)  are  approximated 
by  the  local  mean  and  local  variance  of  the  corrupted  image.  Using  Eqs.  (16),  (17),  and  (15) 
and  applying  static  Kalman  filtering  algorithm,  we  have 


in  which 


Y?  It?  Q;  ■ 

‘,1  2 


An  experimental  example  is  shown  in  Fig.  4,  in  which  the  original  image  is  corrupted 
by  multiplicative  noise  uniformly  distributed  between  0.7  and  1.0,  and  the  estimated  image 
has  been  processed  by  the  algorithm  developed  in  this  section.  Considerable  improvement  is 
shown  in  the  processed  image,  thus  substantiating  the  effectiveness  of  the  local  mean  and 
variance  technique. 


FILTERING  OF  COMBINED  ADDITIVE 
AND  MULTIPLICATIVE  NOISE 

It  is  very  easy  to  extend  the  algorithms  of  previous  sections  to  deal  with  images  cor- 
rupted by  both  additive  and  multiplicative  noise.  A noise-corrupted  image  is  described  by 


^i.y  = ^/.y“i.y  “^/.y 
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(a)  Multiplicative  noise  U(0.7,  1.0)  (b)  Multiplicative  noise  removed 


Figure  4 


in  which  the  statistical  characteristics  are  the  same  as  those  given  in  the  previous  two  sec- 
tions, Assume  that  Uj  j amd  Wj  j are  independent  white  noises.  This  independence  assump- 
tion can  be  removed,  but  the  result  is  a more  complicated  formulation.  Following  the  idea 
of  an  optimal  linear  approximation  of  the  previous  section,  we  have 

The  formulas  for  the  a priori  meam  and  variance  of  Xj  j of  the  previous  section  are  modified 
to  read 

'Xj  j = (Zfj  - J (21) 


and 


9,„,  - 


”2  * “f.J 


‘■J 


(J 


2 
1 ■ 


The  filtering  algorithm  is 


in  which 


^i,>"2  '^1 


(22) 


(23) 
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(a)  Image  corrupted  by  additive  and  (b)  Restored 

multiplicative  noise 


Figure  5 


Figure  5a  shows  the  image  corrupted  by  an  additive  noise  uniformly  distributed  be- 
tween grey  levels  -20  and  -t-20  and  also  a multiplicative  noise  uniformly  distributed  between 
multiplicative  factors  0.7  and  1.0.  The  processed  image,  Figure  5b,  shows  a very  significant 
improvement  over  the  original  image. 


REMARKS  AND  CONCLUSIONS 

The  principal  computational  load  of  the  developed  algorithms  is  in  the  calculation  of 
the  local  means  and  variances,  especially  the  latter.  To  lighten  this  burden,  Wallis  [8]  pro- 
posed a fast  algorithm  in  which  the  im-age  is  partitioned  into  square  subregions  over  which 
the  local  mean  and  variance  are  computed.  Then  the  local  mean  and  variance  of  a particular 
pixel  are  approximated  by  the  use  of  two-dimensional  interpolation  formulas.  Results,  as 
reported  by  Wallis  [8] , are  quite  favorable.  It  is  believed  that  Wallis’  approach  will  yield  an 
equally  impressive  improvement  when  applied  to  our  contrast-enhancement  and  noise- 
filtering algorithms. 

In  conclusion,  image  processing  algorithms  presented  herein,  based  on  considerations 
of  the  local  image  statistics,  have  a structure  that  makes  them  naturally  suitable  for  parallel 
processing.  Since  the  latter  approach  offers  great  computationed  economy,  real-  or  near  real- 
time processing  can  be  achieved.  Future  research  in  this  area  is  to  extend  the  method  to 
image  restoration  of  motion  blur  and  other  degradations  characterized  by  local  correlations 
around  pixels. 
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